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Abstract 

Using both fractional derivatives, defined in the Riemann-Liouville 
and Caputo senses, and classical derivatives of the integer order we ex- 
amine different numerical approaches to ordinary differential equations. 
Generally we formulate some algorithms where four discrete forms of the 
Caputo derivative and three different numerical techniques of solving or- 
dinary differential equations are proposed. We then illustrate how to 
introduce classical initial conditions into equations where the Riemann- 
Liouville derivative is included. 

1 Introduction 

In the past, fractional calculus was applied only from a mathematical point of 
view. The fundamental work was done in ^20 } [2TI f23j 125). At present fractional 
calculus is extremely popular due to a rapid expansion in the field of practical 
applications. Such applications have been used in physics and mechanics [T5tl28). 
finance ^9il26j, hydrology |3J127j and many other disciplines. 

Ordinary differential equations including a mixture of integer and fractional 
derivatives are a natural extension of integer-order differential equations and 
give a novel approach to mathematical modeling many processes in nature. The 
solution of a equation strongly depends on form of the equation and is still 
considered by many authors. It should be noted that an analytical approach is 
Hmited to the Hnear form of equations and includes special functions such as Fox 
and Wright functions O US] or the Mittag-LefHer function [25]. This greatly 
limits practical implementations, i.e. sometimes it is very difficult to illustrate 
the solution in one simple chart. On the other hand, a numerical solution (TJlTl 
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[U [ini \T7\ is an alternative approach to analytical one. However, this approach 
has many disadvantages, i.e. the introduction of the initial conditions included 
in the Riemann-Liouville derivative [l^, the unreasonable assumption that a 
method appHed to a single term equation is proper for solving a multi-term 
equation [HI [21] etc. Against this background Ford [5] noticed that there can 
be a considerable gap between methods that perform well in theory and those 
whose implementations are effective. 

In this paper we try to propose a numerical approach which will be more 
convenient in practical applications. We will give a numerical procedure for how 
to introduce classical initial conditions into an equation where the Riemann- 
Liouville derivative is included. Here we will focus on such types of equation 
as 

/ {x, y (x) , D'yix), D^yix), D^^y{x),..., D"-y{x)) = (1) 

where y(x) is the solution obtained for the class of continuous functions, 
D^y{x), . . . , D^y{x) are derivatives of the integer order, 
D°'^y{x), . . . , D°''^y{x) are derivatives of the fractional order and 
ai, . . . ,am S R are real orders of a fractional derivative. We assume that 
the fractional derivative is defined as the left-side Caputo derivative [2| 

1 ,,(") 

gD"y (x) = — ^ / ^ -^^dT for x > Xq (2) 

Xo 

and the left-side Riemman-Liouville derivative [21] 

.,,D'^y {x) = ^ / ^^^}_,,+, dT for x > xo (3) 

r (n - a) dx" J [x - t) ^ 

Xo 

In above formulae, the notation n — [a] -I- 1 where [• ] is an integer part of a real 
number. Moreover, we introduce a definition of the left-side Riemann-Liouville 
fractional integral [2T| as 

X 

xa 

which will be used in our further calculations. Note that /3 (/3 > 0) is the real 
order of Eqn. (01). On the base of theory [21] we use an expression 

^„Dj2/(x)=,„/r"P"y(x)) (5) 

which shows a relationship between the Caputo derivative ^ and the Riemann- 
Liouville integral ([3]). With regard to papers [Tl|4l[9] in which numerical methods 
are used in the solution of fractional differential equations, the authors mostly 
use the Caputo derivative. However, there is a small number of papers [lO] where 
the authors use the Riemann-Liouville derivative. This small number of papers 
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were confronted by the problem of how to introduce classical initial conditions 
in the Riemann-Liouville derivative in order to obtain a solution for a class of 
continuous functions. In this paper, we propose a way to avoid this problem. 
To be more precise, in every equation where the Riemann-Liouville derivative 
occurs we will change it for the Caputo one. Following this we will discretize 
only the Caputo derivative, except for one case where the real number of the 
Riemann-Liouville derivative dominates in the equation. In this case we propose 
xoD"y{x) = D^xolx'^vi^) ^iid then discretize the left-side Riemann-Liouville 
integral xoIx~°'- 

2 Statement of the problem and its solution 

With regard to Eqn. H]) we limit our considerations to the equation which has 
the following form 

DPy {x) + X< =0 (6) 

where p denotes an integer number being the derivatives order, 
a G (0, 1) is the order of the fractional derivative and A is an arbitrary real 
number. This simple form of the equation allows us to show how our methods 
work properly in comparison to analytical solutions. Note that Eqn. ^ is the 
homogeneous ordinary differential equation with a mixture of derivatives. The 
function y{x) being the solution of this equation, strongly belongs to the class 
of continuous functions. On the basis of our previous results [17] we rewrite 
Eqn. ^ in an expHcit form. Consequently we obtain the three following types 
of equation: 

• p > n for p = 2, a e (0, 1), = 1 

D'y{x) + X'^^D^y{x)^0 (7) 



D^yix) + \x,D^y{x)^0 (8) 

It can be seen in above equations that the integer order of classical deriva- 
tive dominates over the fractional one. 

• p = n for p — I, a d (0, 1), n = 1 

D'y{x) + XZD^y{x)^0 (9) 

D'y{x) + \x„D^y{x)^0 (10) 

In this case we have equal integer orders for the classical and fractional 
derivative. 
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• p < n for p = 0, a £ (0, 1), n = 1 



^^D:y{x)+Xy{x) = (11) 

,,D^y{x) + Xy{x)^0 (12) 

The last case shows fractional ordinary differential equations which are 
well known in the literature. It may observe that the fractional order of 
the equation dominates over the integer one. 

2.1 Analytical solutions 

To compare our direct numerical results we are obligated to solve the above 
system of equations in an analytical way. In this solution we will use a general 
idea which transforms the Riemann-Liouville derivative to the Caputo one [23] . 
Thus we have 

n— 1 / \i~a 

,„D^y (x) = J2 rU_Z + l) ^'y (^o) + ^oD-y (^) (13) 

On the base of |2Q] we also use the Laplace transform. Following that the 
transform of the derivative of the integer order m G iV is 

m— 1 

C [D^y (x)] = s^F (s) ~ s^D^-^-'y {xo) (14) 

fe=0 

The Laplace transform of the Caputo derivative is 

n-1 

C [^Dty {x)] = s^F is) - ^ s'^-'^-'D'^y (xo) (15) 

fc=0 

Using (fTSl and (flSl) we calculated the Laplace transform from the Riemann- 
Liouville derivative in the following form 



CU,D^y{x)]^C 



1=0 ^ ' 



= r?-i+i) + is) - "e s'^-^-'D^y {xo) - ^^^^ 

= s°F{s) 

It should be noted that Eqn. (flGl) is limited by initial conditions which are omit- 
ted here. In previous considerations we assumed function y{x) to be continuous. 
Therefore Eqn. I|16p is contrary to the Laplace transform found in the litera- 
ture [20] where initial conditions of non-integer order occur. This arises from 
an assumption that function y{x) is non-continuous. 
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In our analytical solutions we use also two additional transforms as 

£|j/(i)l = F(») (17) 

Using the above transforms in the set of equations (|7l)- l|12p and retransform- 
ing the results we obtain analytical solutions. Including initial conditions 

y{xQ)^yo, D^y{xa)^y'Q (19) 

the analytical solution to Eqn. ^ is 

y (x) ^ yo + y'a{x~ xq) E2-a,2 (-A [x - xq)^""^ 

D^y [x) = y', E,^^,, (-A {x - xq)'-") ^^^^ 
+ y^, (2 - a) [x - xo) Ei'\, (-A (x - xq)''") 

where i?a,/3 (—Ax") denotes the Mittag-LefHer function [25] which is defined as 

It should be noted that E^^^p (—Ax") occurs in solution l(20|). which is the first 
derivative of the Mittag-Leffler function l(2T1) and is defined as 

E^'^ (-Ax") - V (22) 

i—l ^ ' 

Solving Eqn. ([8]), where initial conditions (|22l) are included, we have 



y (x) = 2/0-B2-a,l (-A (x - Xo)^ 
+ ?/0 (x - Xo) i?2-Q,2 (^-A (x - Xo)^ 



)iy (x) = (2 - a) yoEiXi ("A (x - xo)'"") ^^^^ 

+ J/o E2-a,2 (-A (x - Xo)^"") 

+ (2 - a) y^, (x - xo) E^^, (-A (x - xo)'"") 
Eqn. ([9]) with the initial condition y(xo) = yo has the following solution 

y (x) = yo (24) 

Eqn. (fTO| with the initial condition y(xo) = yo has the following solution 

y (x) = yo£^i-a,i (^-A (x - xo)^~"^ (25) 



Eqn. (fTTI) with the initial condition y{xo) — yo has the following solution 

y{x) = yQEa.i{~X{x-xo)°') (26) 

However Eqn. lfT2|) with the initial condition y{xo) = yo has a trivial solution 
in the following form 

y{x)^0 (27) 

The above solution, arises from our assumption that we only consider a class of 
continuous functions y{x). In our next considerations we neglect this solution. It 
can be seen in literature [9l [10] that the authors solve a similar type of equation 
to Eqn. (fT2|) . where an exciting function f{x) is added on the right side of the 
equation. 



Figure [7] 



2.2 Discrete forms of the Caputo derivative 

There are many propositions on papers [II [131 E] i how to discretize fractional 
operators. Basically several discrete forms are employed to take into account 
the different form of the function included in the fractional derivative. In this 
subsection we would like to propose a general procedure for how to discretize the 
function. Let us consider an independent value x which occurs on a length of 
calculations (xo, xn), where xq and xn are the beginning and end of the range 
respectively. We introduce a homogeneous grid xq < xi < . . . < xn ■ Fig. [T] 
shows four discrete forms of a derivative D"^y{x) = B which is included in the 
Caputo derivative ([2]). 

Taking into account the above discrete forms of the integer derivative we 
propose the following discrete schemes of the Caputo derivative ^ as: 

• the left-side form (case-I) 

(28) 



N 



r(r 



k=l 



(xn - Xk~l) 



{xN - Xk) 



where x £ {xk-i,Xk) and D'^y {xk-i 
• the right-side form (case-II) 



Bk^ 



N 



r(ri- 



k=l 



J2 Bk {xn - Xk-iY^ " - (xn - XkY 



(29) 



where x e {xk-i,Xk) and D"-y{xk) = Bk-, 

• the middle-side form (case-Ill) 

c 



N 



~ T(n-a+l) 51] 



k=l 



{xn ~ Xk-lY' 



where x S {xk-i,Xk) and 



2 



- [xn - XkY' 

Bk+Bk-i 



(30) 
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• the linear form taken from [18] (case-IV) 



where 



JV f 

1 V 1^4^ 

— q) Z-^ I n+1- 

n — Q 



/ \n+l— a / 

(XN - Xk) [XN 



{xn - Xk^lY 



Xk~l 
- [XN - Xk) 



\ n+1 — CK 



(31) 



D"y{xk)-D"y(xk-i) 



Bk ^ D'^yixk) - AkXk. 



We try to use the above forms in numerical schemes to solve ordinary differ- 
ential equations. It should be noted that Eqn. (|28ll (case I) is useful in the ex- 
plicit scheme as, for example, the Euler's method [22] • The next discrete forms 
predicted by formulae l(29|) l|30|) and fSTj) can use for any predictor-corrector 
method [22] , 

Additionally, we also propose a discrete form of the 
Riemann-Liouville fractional integral (01). We consider a range 

(xfe-i, Xk) (fc = 1, . . . , N) for k — 1, . . . , iV and we assume a constant value 
of function y{xk) = Bk. In this case the integral operator |(4]) has the following 
form 



xoIxVi^) 

where Bk — 



1 



r(/3 + i) 



N 



Bi {xN - Xo) -f ^ {Bk - Bk-i) {xN - Xk-i f 



k=2 



(32) 



2.3 Numerical methods of solving ordinary differential equa- 
tions 

In this paper, we chose only three numerical methods which are used in the 
literature to solve an initial-value problem for ordinary differential equations. 

The Euler's method [22] is an explicit one-step method. Using this method 
for any ordinary differential equation of the first order we obtain 

yk = Vk-i + hf {xk-i,yk-i) , k^l,...,N (33) 

The Adams method [22] of the fourth order is the most popular method 
in the literature [5j l6j which is used to solve fractional differential equations. 
This is a predictor-corrector method. It should be noted that for the method 
of fourth order we have to determine three beginning values j/i , ?/2 , 2/3 of the 
function. This can be done by using, for example, the Euler's method. The 
Adams method for a differential equation of the first order has the following 
form: 

• the predictor stage 



• the corrector stage 



f 19 5 1 \ 9 

"''yk = yk-i + h i^fk^i - + ^fk-3 j + ^hfk (35) 

where /fc = /(xfc/'^yfc) , fk-j = f {xk^j ,yk-j) , j = l,...,4, 
fc = 4, ...,iV 

The Gear's method jTTj is also a predictor-corrector type method. This 
method has an advantage in cases where only one call of the function is re- 
quired in one step of the calculation. However, it needs higher derivatives in the 
algebraic form at the beginning point xq. This is a disadvantage of the method. 
The Gear's method for ordinary differential equations of the first order is defined 
as 

• in the predictor stage 

5 

P'^I?'2/fc = ^P,_,D^yfc_i (36) 

3=0 

r for I = -5,..., -1 
where P/ = < ^ , „ ^ and z = 0, . . . , 5, fc = 1, . . . , TV 
1^ for I — 0, ... ,5 

• in the corrector stage 

--D^y^ = P^D^yk + ^^c^A'-D^yk (37) 

where A^^D^yk = "'D^yk - P^^D^yk, '"-D^yk = / {xk,yk,P''D^yk) 
for fc = 1,. ..,7V, i = 0, . ..,5 

and coefficients are defined as 

_ _3_ _ 251 — 1^ — 11 r — — J- 

16 '^1 360' 2 18' 4 6' 5 60' 



All the considered methods will be used in the next sections in order to construct 
algorithms. 

2.4 Algorithms 

In this subsection, we propose several algorithms which solve a set of ordinary 
differential equations presented by general formula 

Algorithm 1.1 

At the beginning we consider Eqn. ^ with initial conditions (fT9|) . On the 
base of the Euler's method l(33|) and the left-side discrete form l(28|) (case-I) of 
Caputo derivative we then propose the following algorithm 
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step 1 Prediction of necessary data: initial conditions, the fractional order a G 
(0, 1), the total length of calculations x £ (xq, xn), the step of calculations 
h. 

step 2 Governing calculations: let k = 1, . . . , TV then 



r(2- 



Uk = Vk-i + hD^Uk-i 



(38) 



Algorithm 1.2 

Considering the same differential equation as presented in Algorithm 1.1 we 
present another approach. This is based on the Adams method l(34|) ([SS)) and 
includes the linear-discrete form l(3T|) (case-IV) of the Caputo derivative. Thus 
we have 

step 1 prediction of necessary data: initial conditions, the fractional order a G 
(0, 1), the total length of calculations x G {xq, xn), the step of calculations 
h and additionally xq^xUo — 0. 

step 2 Introductory calculations: three initial values of the function yi,?/2,2/3 are 
calculated by the Euler's method (algorithm 1.1). 

step 3 Governing calculations: let /c = 4, . . . , TV then the predictor stage is 



^'■y/c = yk-i 



'55 ni 



D'vk-i 



59 ni 



D'yk-2 + ^D'vk- 



,24^ WK-i 24 

P^D^yk = D^yk-i 
iC 

24X0 



9_r,i 

24 



D^Vk-4 



h\ f55C T~ta„, 59C net,, i 37 C ria.. 9 C T]a„, \ 

-'T'^ [Mxo^xVk-l - 24xo^xyk-2 + ^^xo-^xVk-S " 24 xo^a; 2/fc-4 j 



(39) 



l~a) 



(2~a)h 



{D^yj-D^yj-i){xk-Xj)+hD^yj 
{l-a)h 



(Xk-Xj) -(Xk-Xj^i) 



/ \l—a / \1— a 

(Xk-Xj^l) ~ (Xk-Xj) 



(40) 
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and the corrector stage is 



'Vk = Vk-i + h {^D^Vk-i - iiD^yk-2 + j^D'Vk-s) 



'iihD^Vk 



-hX {^c^D^y,_, - + l^lp^Vk^^) - iihXgD^yk 

(41) 

It should be noted that additional assumptions, as shown by step 2 and step 3, 
need to be made for correct calculations to be achieved. 

Algorithm 1.3 

Still considering Eqn. ([7]) with initial conditions ifTO]) we can present another 
approach in comparison to previous ones. This uses the Gear's method l|36p 
l(37|) and the middle-side discrete form (|30l) (case-Ill) of Caputo derivative. 



step 1 Prediction of necessary data: initial conditions, the fractional order a G 
(0, 1), the total length of calculations x £ (xq, xn), the step of calculations 
h and additionally D'^yo = D^yo = D'^yo = = D^yo = 0. 

step 2 Governing calculations: let k — 1, . . . , N then Eqn. ((36l) is the predictor 
stage and beginning of the correction stage is 

^' 1 1 

xo^xVk — r(2-Q) 51] 2 ~ ^i-'^) ~ ^^f^ ~ ^j) 

(42) 

and it follows that we apply the correction stage given by Eqn. (|37|) . 



Algorithm 2.1 

Now let us consider Eqn. ^ with initial conditions (fT9l) . It should be 
noted that this type of equation includes the Riemann-Liouville derivative. We 
apply the Euler's method l(33|) and Eqn. lfT3|) which transforms the Riemann- 
Liouville derivative to the Caputo one. It should be remembered that Eqn. ifTS)) 



n-l f^_^ y-a 

includes a set of initial conditions in the general form ^ rd-a+i) (^o)- 

In the case of the lower Hmit x tends to a;o and then the above expression is 
infinite. This means that we cannot use a full numerical approach in order to 
solve the class of ordinary differential equations where the Riemann-Liouville 
derivative is included. However an assumption of homogeneous initial condi- 
tions i.e. D^y [xq) = avoids the problem. Eqn. has a correct analytical 
solution (|23l) and the function y[x) has a finite value at the beginning point xq 
which is contrary to previous considerations. With regard to our assumption 
that function y{x) belongs to the class of continuous functions we are obHgated 
to improve the factor (x — x^f " included in Eqn. (fT3|) in order to restrict the 
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function continuity. Therefore we put {x* — xqY " into Eqn. ifTSj) instead of 

(x — xqY " where x* — x + c and c — {T {1 — a)) " . Point c is where the shift 
from initial conditions to homogeneous initial conditions occurs. Our assump- 
tion allows for the correct behaviour of factor {x* — a^o)* " — {x + c — xq)'"" 
in Eqn. ifTSj) in the class of continuous functions especially when x tends to Xq. 
Taking into account the above considerations we have the following algorithm 

step 1 prediction of necessary data: initial conditions, the fractional order 
a £ (0, 1), the total length of calculations x G (xq^xn), the step of calcu- 
lations h. 

step 2 Governing calculations: let k = 1, . . . ,N then 

+ r(2-a) E D^Vk-l {Xk - Xj-i) " ~{xk- Xj) 



(43) 



Vk = J/fc-i + hD^yk^i 

where x^_-^ = Xk-i + he. 
Algorithm 2.2 

The Euler's method has been used in algorithm 2.1. Considering Eqn. ^ 
with initial conditions lfT9|) and using the Adams method l(34|) with the linear- 
discrete form (|3T|) (case-IV) of Caputo derivative we obtain 

step 1 prediction of necessary data: initial conditions, the fractional order a G 
(0, 1), the total length of calculations x £ (a;o, xn), the step of calculations 

h, 2^D<^yo = and c = (r (1 - a))""". 

step 2 Introductory calculations: three initial values of the function yi,y2,y3 are 
calculated by the Euler's method (algorithm 2.1). 

step 3 Governing calculations: let A: = 4, . . . , iV then the predictor stage is 

P^yk = yk-i + h (fi^iyfc-i - ^D^yk-2 + ^D^Vk^^ - ^D^yk-i) 

P^D^yk = D^yk^i 
-h\ {^x^D'^yk-i - ^xoD"yk^2 + ^xoDxUk-s - ^xoD^Vk-A) 

(44) 
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lo^xlJk — r(l-Q) 



-yo 



r( 



{2-a)h 



{Xk-Xj) -(Xk-Xj^i) 



+ 



(D^yj-D\j^i){xk-Xj)+hD^yj 
{l-a)h 



{Xk -Xj-i) 



1-Q 



(Xfc - Xj) 



1-a 



(45) 



and the corrector stage is 



^hXxfjD'^yk 



(46) 



Algorithm 2.3 

We present another algorithm on the base of Algorithms 2.1 and 2.2. We still 
consider Eqn. ^ with the initial conditions lfT9|) and then we use the Gear's 
method (|36|) l(37l) with the middle-side discrete form (|3Q|) (case-Ill) of Caputo 
derivative. We have 

step 1 prediction of necessary data: initial conditions, the fractional order a g 
(0, 1), the total length of calculations x £ (xq, a; at), the step of calculations 

/i, c = (r (1 - a))-"'' and D^yo = D^yo = = D^yo = D^yo = 0. 



step 2 Governing calculations: let k = 1,. . . ,N then Eqn. (|36|) is the predictor 
stage and beginning of the correction stage is 



D^yk+D^yt-i 

-{2-a) ^ 2 



^^xVk — r(2-Q) S 



/ \l — a f \\ — OL 

(Xfe-Xj-i) -{Xk-Xj) 



(47) 



and it follows that we apply the correction stage given by Eqn. (|37|) . 



Algorithm 3.1 

Next we consider another class of ordinary differential equations given by 
formula ([9]) with the initial condition yix^) = yo. Using the Euler's method l|33p 
and the left-side discrete form ((28l) (case-I) of the Caputo derivative we then 
propose the following algorithm 

step 1 prediction of necessary data: initial conditions, the fractional order a € 
(0, 1), the total length of calculations x G {xo,xn), the step of calculations 
h and D^yo — 0. 
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step 2 Governing calculations: let k = 1, . . . ,N then 



(48) 



k 

I 

Vk = yfc-i - hX^^D^yk-i 
D'vk = -X^D^Vk 

To simplify our considerations we neglect other algorithms which can be 
applied to solve Eqn. ([9]). On the basis of previous explanations of predictor- 
corrector methods one can construct interesting procedures self-reHantly. 

Algorithm 4.1 

Using the Euler's method ((33l) and the left-side discrete form l|28p (case-I) of 
the Caputo derivative we construct an algorithm which solves Eqn. lfTO|) with 
the initial condition y (xo) — j/o- Thus we obtain 

step 1 prediction of necessary data: initial conditions, the fractional order a G 
(0, 1), the total length of calculations x £ {xo,xn), the step of calculations 
h, D^yo = and c = (r (1 - a))"""' . 

step 2 Governing calculations: let k = 1, . . . ,N then 

xa^xVk - r(l-a) 

+ r(2-a) E D^Vk^l (xk - Xj-i) " - {Xk~ Xj) 



(49) 

yk = yk-i - hX^^„D^yk-i 



D^yk = -K.D'^yk 

where x*f^_i ~ Xk-i + he. 
We also neglect other constructions of algorithms. 

Algorithm 5.1 

Using the Euler's method ((33l) we solve numerically the last equation 
taking into account the initial condition y (xq) = yo- Follow the discrete form 
of the Riemann-Liouville fractional integral (|32l) we have 

step 1 prediction of necessary data: initial conditions, the fractional order a € 
(0, 1), the total length of calculations x £ {xo,xn), the step of calculations 
h. 
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step 2 Governing calculations: let k = 1, . . . ,N then 

- I] -Sj-i (2;fc-a;j-2) (xk-Xj^i) " ^ (50) 
yfc = Vk-i + hBk 

In summary we propose numerical schemes suitable for all possible cases 
of Eqn. We then generally apply two approaches. One is connected with 
the explicit method and the second is based on the predictor-corrector scheme. 
It should be noted that we illustrate all the considered methods (the Euler's 
method, the Adams method and the Gear's method) which are used in Eqns |[7]) 
([HI. However we limit this illustration in Eqn. f9|- lfTT|) only to the Euler's 
method. In this case we omit the predictor-corrector methods because there 
are problems with applications. Using this approach we found that the Adams 
method requires known values of the first derivative in the algebraic form. How- 
ever the Gear's method requires the algebraic form of the second derivative 
which is used in calculations of the correction factor. 



3 Results 

The algorithms presented in the previous section allow us to compare numerical 
and analytical results respectively. First we compare the results obtained from 
four discrete forms of the Caputo derivative. Thus we assume a function 

y i^) = (51) 

where the Caputo derivative is 

Table[T]shows analytical values of the function l|52p at assumed points. The other TafeZefT] 
rows of this table present the difference between the numerical and analytical 
results. Analysing this table we can see that the liner-discrete form of the 
Caputo derivative ((3T|) (case-IV) gives the best results. However, this discrete 
form is very complex. This is a disadvantage of case-IV. It should be noted that 
the middle-side form (|30l) (case-Ill) also gives quite small errors. The next two 
discrete forms l(28|) and l(29l) generate large errors in comparison to the previous 
cases. 

We try to estimate how several discrete forms of the Caputo derivative work 
in numerical schemes used to solve ordinary differential equations. It should be 
remembered that, as presented in subsection 2.4, there are many constructions 
of numerical approaches which are dependent on the form of the ordinary dif- 
ferential equation, the discrete form of the Caputo derivative assumed (cases I 
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to IV) and, finally, on the numerical scheme (the explicit or predictor-correct 
schemes) used. Therefore, a general question arises: how to efficiently compare 
all the possible constructions of numerical schemes in order to comment on their 
practical application? First of all, we started our comparison between predictor- 
corrector methods. We can apply three forms of the Caputo derivative (cases II 
to IV) and two numerical schemes (the Adams and Gear's schemes) as presented 
in the previous section. Let us start with Eqn. ^ where initial conditions 

y(0) = 0, D'y {0) = 1 (53) 

are assumed. Note that Eqn. ([7]) has an analytical solution ((2Ql) . 
Table [H Table [2] presents a direct comparison of the two predictor-corrector meth- 

ods where three discrete forms of the Caputo derivative for each method are 
included. It should be noted that we obtained many sets of data. Brief anal- 
ysis presents the best results in the linear-discrete form (case-IV) fSTj) of the 
Caputo derivative for both the Gear's l(36|) l(37l) and Adams fM)) (|35l) methods. 
Satisfactory results can be also achieved by the application of the middle-side 
form (case-Ill) (|30l) obviously used for both methods. Taking into account the 
complexity of the forms, we chose the middle-side form ((30|) for the next calcu- 
lations using the predictor-corrector method. We also noticed that the Gear's 
method l(36|) (|37|) generates smaller error values than the Adams method l(34l) 
l(35|l. 

y(0) = -l, D'y {0)^1 (54) 

Now we consider all the numerical methods used to solve the same equation |[7]) 
Figure [H with initial conditions l(54|) . Fig. [2] illustrates the analytical solution taking into 

account the influence of parameter a. Notice that for such a solution we used 
Table[3l algorithms 1.1, 1.2 and 1.3 respectively. Table [3] presents the analytical values 

given by formula l(20l) and errors given by several algorithms used. It can be 
observed that the Gear's method (algorithm 1.3) generates the smallest errors 
in comparison to other methods used. However, the Euler's method (algorithm 
1.1) performs better than the Adams method (algorithm 1.2). As we noted 
in the previous section, predictor-corrector methods have some disadvantages 
which are revealed in the calculations of higher derivatives for the Gear's method 
or necessary values of the function at the initial three points for the Adams 
method. Moreover, a discrete form of the Caputo derivative has less of an 
influence on error calculations than a numerical method used in solving an 
ordinary differential equation. 

Next we consider all the numerical methods used to solve of Eqn. ^ with 
initial conditions 

y(0) = l, D'y {0) = 1 (55) 

This equation has an analytical solution given by formula l|23p It should be 

noted that for a such solution we used algorithms 2.1, 2.2 and 2.3 respectively. 

Fig. [3] shows the behavior of the analytical solution over an independent value x Figure[M 

for different values of the parameter a. Table [4] presents analytical values given Taft/ej^] 

by formula l(23|) and errors given by several algorithms used. This numerical 
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solution confirms our previous considerations that the Euler's and Gear's meth- 
ods show the smallest errors. Taking into account both the errors generated by 
the numerical method and the complexity of numerical schemes we chose the 
Euler's method as the method to use in the next calculations. 

Assuming the initial condition y (0) = 2 for Eqn. ^ we obtain an analytical 
solution in the form 

y{x)^2 (56) 

Fig. [4] shows direct comparison between analytical solution l(56|) and the Figure\4\ 
Euler's method (algorithm 3.1) for Eqn. ([9]). In this case we have not printed 
a table because only a very small difference between the analytical and numerical 
solutions can be observed. 

Next we consider Eqn. (fTO| with the initial condition y(0) = 1. This equation 
has an analytical solution given by formula (|25|) which is presented in Fig.[5l We Figure[M 
then apply three numerical schemes in order to compare data with the analytical 
results. Table [5] shows the analytical values given by formula (|25|) and errors Table[^ 
given by the Euler's method (algorithm 4.1). It may be observed that errors 
generated by the Euler's method are small. This confirms that the Euler's 
method is an adequate method for solving ordinary differential equations with 
the a mixture of derivatives. 

The last case concerns the fractional differential equation given by formula 
ifTTj) . We assume the initial condition y{0) — 1. This equation has an analytical 
solution presented by formula ((26l) . We use the Euler's method given by algo- 
rithm 5.1 for a numerical solution. Fig. [6] presents the analytical solution l(26l) Figure[d[ 
over an independent value x for different values of parameter a. Table [6] shows Table\di 
analytical values and errors generated by the Euler's method (algorithm 5.1). It 
can be observed that the Euler's method (algorithm 5.1) generates small error 
values. Summarising our results we proved that the Euler's method is suitable 
for the numerical solution of ordinary differential equations having a mixture 
of derivatives. This method has the following advantages: simpHcity, a small 
difference between analytical and numerical values depending on the step of 
calculations, stable error values over all the considered length of calculations. 



4 Conclusions 

In this study, we proposed numerical algorithms to solve ordinary differential 
equations where the a mixture of fractional- and integer-order derivatives oc- 
curs. We used three known numerical techniques, the Euler's, Adams and Gear's 
methods, to solve such equations. Taking into account the equation order, we 
divided differential equations into three classes; where the integer order dom- 
inated over an integer number calculated from the fractional order, where the 
integer order and number were the same and where the integer number domi- 
nated over the integer order. In the considered equations we distinguished two 
types of fractional operators: the left-side Riemann-Liouville and left-side Ca- 
puto derivatives. Using a known transition rule between the derivatives, in a 
differential equation where the Riemman-Liouville operator occurs, we changed 
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the Riemann-Liouville derivative to the Caputo derivative. Next we proposed 
four discrete forms of the Caputo derivative. On the basis of the previous classifi- 
cation of ordinary differential equations we illustrated the proper algorithms. It 
should be noted that our algorithms are valid in the class of continuous functions. 
This assumption allows us to solve the problem of how to include classical ini- 
tial conditions into ordinary differential equations where the Riemann-Liouville 
fractional derivative occurs. 

Using direct comparison between the analytical and numerical data we ob- 
tained satisfactory results. Deeper analysis shows that the 
predictor-corrector methods (the Adams and Gear's methods) require the linear- 
discrete form (|3T|) of the Caputo derivative in order to refiect the analytical data 
more precisely. Satisfactory results can also be achieved by the application of 
the middle-side form (|30| obviously used for both the methods. Comparing the 
results obtained by the above methods we observed that the Gear's method gives 
better results than the Adams method. We also compare the results obtained by 
the analytical solution, the predictor-corrector methods and the Euler's method 
which we first successfully applied to solve ordinary differential equations in- 
cluding mixture of derivatives. We can say that the Euler's and Gear's methods 
show the smallest errors. This may be unexpected because the Euler's method 
is a method of the order O (h) and additionally it includes the left-side discrete 
form (|28|) of the Caputo derivative. On the other hand the predictor-corrector 
methods require additional circumstances, for example the Adams method of 
fourth order needs the four initial values of the function to be determined and 
the Gear's method needs higher derivatives in algebraic form at the starting 
point xq. There are disadvantages of predictor-corrector methods. Against this 
background we observed that a discrete form of the Caputo derivative and the 
method order have less of an infiuence on error calculations than the disadvan- 
tages of a numerical method in solving an ordinary differential equation. Taking 
into account both the errors generated by the numerical method and the com- 
plexity of numerical schemes we chose the Euler's method as a suitable method 
to use for practical calculations. 
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List of captions for illustrations 



1 Discrete forms of an integer derivative D'^y{x) — B for the range {xk-i, Xk) 
for fc= 1,...,7V: 

a) left-side 

b) right-side 

c) middle-side 

d) linear 

2 Analytical solution of Eqn. ([7]) with initial conditions l(54l) over an inde- 
pendent value X for different values of a. 

3 Analytical solution of Eqn. ([8]) with initial conditions l|55p over an inde- 
pendent value X for different values of a. 

4 Comparison of analytical l(56|) and numerical (algorithm 3.1) for Eqn. |[9]) 
with the initial condition y{0) = 2. 

5 Analytical solution of Eqn. ifTO|) with the initial condition 
2/(0) = 1 over an independent value x for different values of a. 

6 Analytical solution of Eqn. ifTTj) with the initial condition 
2/(0) = 1 over an independent value x for different values of a. 
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List of captions for tables 



Table [U Analytical values of the Caputo derivative (|52l) and errors generated by 
discrete forms of the Caputo derivative ^ ^ §0^ ^ for h = 0.01. 

Table E] Analytical solution of Eqn. ([7]) with initial conditions l(53l) and errors 
generated by the Gear's and Adams methods for h — 0.01. 

Table [3] Analytical results of Eqn. ^ with initial conditions (|54l) and errors gen- 
erated by the Euler's (algorithm 1.1), Adams (algorithm 1.2) and Gear's 
(algorithm 1.3) methods for h = 0.01. 

Table [4] Analytical results of Eqn. |[8]) with initial conditions (|55l) and errors gen- 
erated by the Euler's (algorithm 2.1), Adams (algorithm 2.2) and Gear's 
(algorithm 2.3) methods for h = 0.01. 

Tabled] Analytical results of Eqn. ifTOl) with the initial condition 

y{0) = 1 and errors generated by the Euler's method (algorithm 4.1) 
for h = 0.01. 

Table [6] Analytical results of Eqn. ifTTj) with the initial condition 

y(0) — 1 and errors generated by the Euler's method (algorithm 5.1) 
for h = 0.01. 
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Table 1: 





2/(1) 


y(4) 


2/(6) 


2/(8) 


2/(10) 


a = 0.1 


analytical 


1.0944780 


15.2447755 


32.9377877 


56.8955141 


86.9374806 


errors generated by discrete forms of the Caputo derivative 


case-I 


1.04e-2 


3.62e-2 


5.22e-2 


6.76e-2 


8.26e-2 


case-II 


1.04e-2 


3.62e-2 


5.21e-2 


6.75e-2 


8.26e-2 


case-Ill 


1.77e-5 


1.98e-5 


2.02e-5 


2.07e-5 


2.10e-5 


case-IV 








2.00e-8 





4.00e-8 


a = 0.5 


analytical 


1.5045056 


12.0360444 


22.1116256 


34.0430746 


47.5766431 


errors genei 


■ated by discrete forms of the Caputo derivative 




case-I 


1.17e-2 


2.30e-2 


2.81e-2 


3.24e-2 


3.61e-2 


case-II 


1.08e-2 


2.21e-2 


2.72e-2 


3.14e-2 


3.52e-2 


casc-III 


4.60C-4 


4.64e-4 


4.65e-4 


4.660-4 


4.66e-4 


case-IV 











2.00e-8 


l.OOe-8 


a = 0.9 


analytical 


1.9111582 


8.7813771 


13.7171223 


18.8232939 


24.0600562 


errors generated by discrete forms of the Caputo derivative 


case-I 


1.60e-2 


1.76e-2 


1.81e-2 


1.8oe-2 


1.88e-2 


caso-II 


4.980-3 


C..540-3 


7.04O-3 


7.410-3 


7.70O-3 


case-Ill 


o.53e-3 


5.53e-3 


5.53e-3 


5.53e-3 


5.o3e-3 


case-IV 











l.OOe-8 


l.OOe-8 
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Table 2: 



2/(1) 2/(4) 2/(6) 


2/(8) 


2/(10) 


a = 0.1 


analytical 0.8226218 -0.4960927 -0.2222706 


0.5587627 


-0.1911696 


errors generated by the Gear's method 


case-II 2.71e-4 1.39e-2 4.21e-3 


5.77e-3 


1.39e-2 


case-Ill 1.40e-5 1.24e-5 4.23e-6 


1.47e-5 


7.66e-6 


case-IV 1.06e-5 1.7oe-5 3.08e-5 


5.57e-6 


3.18e-5 


errors generated by the Adams method 


case-II 8.55e-5 1.40e-2 4.27e-3 


5.89e-3 


1.40e-2 


case-Ill 1.71e-4 1.05e-4 5.35c-5 


l.llc-4 


3.97e-5 


case-IV 1.70e-4 1.14e-4 3.67e-5 


1.09e-4 


5.85e-5 


a = 0.5 


analytical 0.7374822 0.3129-516 0.16239.5.5 


0.2017797 


0.18C7275 


errors generated by the Gear's method 


case-II 7.85e-4 7.84e-3 5.76e-3 


4.45e-3 


5.04e-3 


case-Ill 1.26e-4 1.02e-4 1.96e-5 


2.53e-5 


8.40e-6 


casc-IV 1.67C-4 7.25c-5 4.27e-5 


4.71e-5 


4.29e-5 


error generated by the Adams method 


case-II 6.92e-5 7.54e-3 5.62e-3 


4.27e-3 


4.87e-3 


case-Ill 5.90e-4 3.97e-4 1.22e-4 


1.56C-4 


1.75C-4 


case-IV 5.40e-4 2.22e-4 1.06e-4 


1.34e-4 


1.23e-4 


a = 0.9 


analytical 0.6512921 0.8643086 0.8128471 


0.7787088 


0.7563200 


errors generated by the Gear's method 


case-II 2.29e-3 4.97e-3 5.13e-3 


5.09e-3 


5.06e-3 


case-Ill 1.09e-3 5.18e-5 1.20e-4 


7.50e-5 


4.63e-5 


case-IV 1.80C-3 2.39c-3 2.25e-3 


2.15e-3 


2.09e-3 


errors generated by the Adams method 


case-II 8.76e-5 1.78e-3 2.13e-3 


2.23e-3 


2.28e-3 


casc-III 1.25e-3 3.19e-3 3.06e-3 


2.89e-3 


2.78e-3 


case-IV 5.65e-4 7.89e-4 7.37e-4 


7.02e-4 


6.79e-4 
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Table 3: 







2/(4) 


2/(6) 


2/(8) 


2/(10) 


a = 0.1 


analytical 


-0.1773782 


-1.4960927 


-1.2222706 


-0.4412373 


-1.1911696 


Euler 


1.79e-5 


2.68e-5 


2.92e-5 


4.23e-6 


3.85e-5 


Gear 


1.40e-5 


1.24e-5 


4.18e-6 


1.47e-5 


7.60e-6 


Adams 


1.70e-4 


1.14e-4 


3.67e-5 


1.09e-4 


5.85e-5 


a = 0.5 


analytical 


-0.2625177 


-0.6870484 


-0.8376044 


-0.7982203 


-0.8132725 


Eulcr 


1.26e-4 


1.04e-4 


2.41e-5 


2.51e-5 


9.32e-6 


Gear 


1.26e-4 


1.02e-4 


1.96e-5 


2.53e-5 


8.40e-6 


Adams 


5.40e-4 


2.22e-4 


1.06e-4 


1.34e-4 


1.23e-4 


a = 0.9 


analytical 


-0.3487079 


-0.1356914 


-0.1871529 


-0.2212912 


-0.2436800 


Euler 


1.08e-3 


5.29e-5 


1.20e-4 


7.49e-5 


4.62e-5 


Gear 


1.06e-3 


5.18e-5 


1.20e-4 


7.50e-5 


4.63e-5 


Adams 


5.65e-4 


7.89e-4 


7.37e-4 


7.02e-4 


6.79e-4 



Table 4: 





2/(1) 


2/(4) 


2/(6) 


2/(8) 


2/(10) 


a = 0.1 


analytical 


1.3290813 


-1.0029826 


0.3872104 


0.4928008 


-0.5873420 


Euler 


4.83e-3 


3.40e-3 


8.51e-4 


3.06e-3 


1.49e-3 


Gear 


4.83e-3 


3.41e-3 


8.92e-4 


3.10e-3 


1.49e-3 


Adams 


1.06e-2 


6.39e-3 


2.50e-3 


7.28e-3 


2.39e-3 


a = 0.5 


analytical 


1.1341116 


0.1100800 


0.1748923 


0.2090891 


0.1714270 


Euler 


7.42e-3 


8.06e-5 


3.56e-4 


8.41e-4 


6.45e-4 


Gear 


7.42e-3 


7.54e-5 


3.53e-4 


8.43e-4 


6.45e-4 


Adams 


2.27e-2 


7.88e-3 


4.71e-3 


6.06e-3 


5.61e-3 


a = 0.9 


analytical 


1.0146791 


0.8374876 


0.7913672 


0.7652525 




Euler 


5.73e-3 


2.17e-3 


2.01e-3 


1.96C-3 




Gear 


5.74e-3 


2.17e-3 


2.01e-3 


1.96e-3 




Adams 


3.15e-2 


3.69e-2 


3.50e-2 


3.37e-2 
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Table 5: 





y(i) 


2/(2) 


y{3) 


y(4) 


2/(5) 


a = 0.1 


analytical 


0.3760660 


0.1811155 


0.1014866 


0.0644356 


0.0452231 


Euler 


1.36e-3 


9.40e-4 


5.42e-4 


3.03e-4 


1.70e-4 


n = ()..■) 


analytical 


0.4275836 


0.3362040 


0.2873412 


0.2553957 


0.2323262 


Euler 


9.61e-4 


8.60e-4 


7.88e-4 


7.31e-4 


6.85e-4 


a = 0.9 


analytical 


0.4855645 


0.4682030 


0.4580801 


0.4509182 


0.4453768 



Euler 1.23e-3 1.19e-3 1.17e-3 1.16e-3 1.14e-3 



Table 6: 





2/(1) 


2/(2) 


2/(3) 


2/(4) 


2/(5) 


a = 0.1 


analytical 


0.4855645 


0.4682030 


0.4580801 


0.4509182 


0.4453768 


Euler 


1.33e-4 


7.06e-5 


4.85e-5 


3.69e-5 


3.07e-5 


a = 0.5 


analytical 


0.4275836 


0.3362040 


0.2873412 


0.2553957 


0.2323262 


Euler 


8.39e-4 


5.30e-4 


3.73e-4 


2.82e-4 


2.23e-4 


a = 0.9 


analytical 


0.3760660 


0.1811155 


0.1014866 


0.0644356 


0.0452231 


Euler 


1.55e-3 


1.27e-3 


7.99e-4 


4.77e-4 


2.88e-4 
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